www.gusucode.com > A benchmark software for MSPC工具箱matlab程序 > A benchmark software for MSPC/qFactor.m
function [sys,x0,str,ts] = qFactor(t,x,u,flag) %qFactor_S_Function: Calculates q factor from composition z and temperature T switch flag, %%%%%%%%%%%%%%%%%% % Initialization % %%%%%%%%%%%%%%%%%% case 0, [sys,x0,str,ts]=mdlInitializeSizes(); %%%%%%%%%%%%%%% % Derivatives % %%%%%%%%%%%%%%% case 1, sys=mdlDerivatives(t,x,u); %%%%%%%%%%% % Outputs % %%%%%%%%%%% case 3, sys=mdlOutputs(t,x,u); %%%%%%%%%%%%%%%%%%% % Unhandled flags % %%%%%%%%%%%%%%%%%%% case { 2, 4, 9 }, sys = []; %%%%%%%%%%%%%%%%%%%% % Unexpected flags % %%%%%%%%%%%%%%%%%%%% otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end % end csfunc % %============================================================================= % mdlInitializeSizes % Return the sizes, initial conditions, and sample times for the S-function. %============================================================================= % function [sys,x0,str,ts]=mdlInitializeSizes() sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 2; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); % % initialize the initial conditions % x0 = []; % % str is always an empty matrix % str = []; % % initialize the array of sample times % ts = [0 0]; % Specify the block simStateCompliance. The allowed values are: % 'UnknownSimState', < The default setting; warn and assume DefaultSimState % 'DefaultSimState', < Same sim state as a built-in block % 'HasNoSimState', < No sim state % 'DisallowSimState' < Error out when saving or restoring the model sim % state % end mdlInitializeSizes % %============================================================================= % mdlDerivatives % Return the derivatives for the continuous states. %============================================================================= % function sys=mdlDerivatives(t,x,u) sys = []; % end mdlDerivatives % %============================================================================= % mdlOutputs % Return the block outputs. %============================================================================= % function sys=mdlOutputs(t,x,u) %Data %******************************** VLEdata = load('VLEdata.mat'); xA = VLEdata.xA; Tbubble = VLEdata.Tbubble; Tdew = VLEdata.Tdew; M_methanol = VLEdata.M_methanol; M_ethanol = VLEdata.M_ethanol; CpL_methanol = VLEdata.CpL_methanol; CpL_ethanol = VLEdata.CpL_ethanol; CpV_methanol = VLEdata.CpV_methanol; CpV_ethanol = VLEdata.CpV_ethanol; Tn_methanol = VLEdata.Tn_methanol; Tn_ethanol = VLEdata.Tn_ethanol; Tc_methanol = VLEdata.Tc_methanol; Tc_ethanol = VLEdata.Tc_ethanol; Hv_methanol = VLEdata.Hv_methanol; Hv_ethanol = VLEdata.Hv_ethanol; %******************************** %Calculate Tbubble and Tdew from zF value: z = u(1); T = u(2); Tbp = interp1(xA,Tbubble,z); Tdp = interp1(xA,Tdew,z); %Possible cases: % Feed condition q Case %--------------- ----- ------------ % Subcooled liquid >1 (1) TF < Tb % Saturated liquid 1 (2) TF = Tb % Vapor-liquid mixture 0<q<1 (3) Tb < TF < Td % Saturated vapor 0 (4) TF = Td % Superheated vapor <0 (5) TF > Td if T<Tbp; %Subcooled liquid CpL = CpL_methanol*z+CpL_ethanol*(1-z); %MF = M_methanol*z+M_ethanol*(1-z); HL = CpL*(Tbp-T); %Corrección por temperatura: Hv_methanol = Hv_methanol*1000*((Tc_methanol-(Tbp+273.15))/(Tc_methanol-Tn_methanol))^0.38; Hv_ethanol = Hv_ethanol*1000*((Tc_ethanol-(Tbp+273.15))/(Tc_ethanol-Tn_ethanol))^0.38; HG = z*(CpL_methanol*(Tdp-T)+Hv_methanol)+(1-z)*(CpL_ethanol*(Tdp-T)+Hv_ethanol); q = HG/(HG-HL); elseif T==Tbp; %Saturated liquid q = 1; elseif Tbp<T && T<Tdp %Vapor-liquid mixture elseif T==Tdp; %Saturated vapor q = 0; elseif T>Tdp %Superheated vapor CpV = CpV_methanol*z+CpV_ethanol*(1-z); %Corrección por temperatura: Hv_methanol = Hv_methanol*((Tc_methanol-T)/(Tc_methanol-Tn_methanol))^0.38; Hv_ethanol = Hv_ethanol*((Tc_ethanol-T)/(Tc_ethanol-Tn_ethanol))^0.38; Hv = Hv_methanol*z+Hv_ethanol*(1-z); q = -CpV*(T-Tdp)/Hv; end sys = q; % end mdlOutputs